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Abstract: This paper presents a two stage algorithm for real-time estimation of instanta- 
neous tremor parameters from gyroscope recordings. Gyroscopes possess the advantage of 
providing directly joint rotational speed, overcoming the limitations of traditional tremor 
recording based on accelerometers. The proposed algorithm first extracts tremor patterns 
from raw angular data, and afterwards estimates its instantaneous amplitude and frequency. 
Real-time separation of voluntary and tremorous motion relies on their different frequency 
contents, whereas tremor modelling is based on an adaptive LMS algorithm and a Kalman 
filter. Tremor parameters will be employed to drive a neuroprosthesis for tremor suppression 
based on biomechanical loading. 

Keywords: tremor; inertial sensors; MEMS gyroscope; tremor modelling; voluntary move- 
ment estimation; adaptive signal processing; Kalman filter; real-time estimation; neuropros- 
thesis 



1. Introduction 

Tremor is defined as a rhythmic oscillation of a body part [1]. Despite we all have a small component 
of tremor, the so-called physiological tremor, there exist pathologies with very disabling forms of tremor. 
Pathological tremor constitutes the most common movement disorder and is continuously increasing its 
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prevalence with ageing. Although not life threatening, tremor often causes social embarrassment and 
functional disability; in fact, 65% of patients suffering from upper limb tremor report serious difficulties 
in performing their activities of daily living [2]. 

Regarding its aetiology, tremors are caused by the different combinations of four different mech- 
anisms: (1) mechanically induced oscillations, (2) oscillations due to reflexes, (3) oscillations due to 
central neuronal circuits, and (4) oscillations because of disturbed feedback or feedforward loops [1]. 
However, understanding of the origins of tremors is far from complete. 

Pathological tremors are typically classified by position/motor behavior, which reflects under which 
circumstances tremor appears. According to this criterium, tremor falls into three categories: rest, pos- 
tural and kinetic tremor [3]. 

Current strategies on the treatment of tremors are based on drugs, surgery (thalatomy), or deep brain 
stimulation. The last two options are typically employed in patients refractory to drugs [4]. How- 
ever, (1) tremor is not managed satisfactorily in 25 % of patients, (2) drugs used often induce side 
effects, and (3) neurosurgery is associated with hemorrhages and psychiatric manifestations, such as an 
increased suicidal tendency [5]. These reasons make research on new therapeutic options to manage 
tremor mandatory. 

In this context, functional compensation of pathological tremors via biomechanical loading has ap- 
peared as a promising alternative. It is based on the fact that most types of tremors change their char- 
acteristics (amplitude and frequency) when the apparent limb impedance is modified, for example, ap- 
plying a force or adding a mass [6-8]. In particular, it has been clinically demonstrated that increase 
of damping and/or inertia in the upper limb leads to an effective tremor reduction [2,8]. In [9], the 
authors have validated, both functionally and clinically, the concept of tremor suppression by means 
of a wearable robot that applies biomechanical loads in the upper limb. To achieve satisfactory tremor 
suppression, the wearable robot must leave concomitant voluntary movement unmodified, compensating 
for tremor. In this regard, real-time estimation of instantaneous tremor parameters becomes fundamen- 
tal to develop robust orthotic or neuroprosthetic solutions for tremor suppression. The present work is 
carried out in the framework of EU project TREMOR (ICT-2007-224051), which aims at developing 
a soft wearable robot for tremor suppression based on a textile substrate. Tremor suppression will be 
achieved through the application of selective biomechanical loads by means of Functional Electrical 
Stimulation (FES). FES constitutes a means of activating motoneurons or reflex pathways by stimu- 
lating sensory nerve fibers [10]; it therefore allows for employing human muscles as actuators for a 
wearable robot [11]. The algorithm presented here will be employed to derive tremor characteristics, 
in order to generate adequate stimulation patterns that counteract tremor, leaving concomitant voluntary 
movement unaffected. 

Inertial sensors technologies constitute a breakthrough and optimal approach to worn sensors in mo- 
tion analysis [12], and wearable robotics [13]. Advances in microelectromechanical systems (MEMS) 
permitted to develop miniaturized low cost sensors that measure changes in velocity, position and accel- 
eration [13]. Low size and weight allow for wearability, overcoming the limitation of traditional motion 
capture systems, based on optic sensors and markers, to work in reduced areas. The most spread types of 
MEMS sensors are solid state accelerometers, gyroscopes, and magnetometers. MEMS accelerometers, 
and most recently gyroscopes, have been employed to analyze tremorous movements. 
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Other solutions for kinematic recording of tremors found in the literature are instrumented, ad-hoc 
devices, such as the instrumented glove [14], or digitizing tablets [15]. Works that study tremor by 
generating it with an actuated mechanical setup, employ specific sensors technologies such as Hall effect 
transducers [16], or piezoelectric accelerometers [17], to measure it. 

Regarding tremor monitoring with MEMS inertial sensors, accelerometers constitute the most spread 
alternative [18,19], although it is not the most adequate, [20]. In fact, accelerometers measure limb ori- 
entation with respect to gravity, but their measurement is corrupted by voluntary movements, decreasing 
accuracy dramatically. Moreover, accelerometers measure linear acceleration, whereas articular motions 
are rotational about joints. Accelerometer data is typically considered to be compound of three factors: 
linear acceleration, gravity, and white additive noise [21]. No analytic model to distinguish between 
accelerometer data due to acceleration and gravity is available, although the most common approach is 
to model acceleration changes as a low pass filter [21,22], or to separate both components by means of 
a FIR low pass filter [23]. On the contrary, solid state gyroscopes provide a direct measurement of a 
rotational movement, uninfluenced by gravity [24], but affected by a characteristic low frequency bias. 
This bias is traditionally thought of being caused by the effects of temperature and mechanical wear [21], 
although our experiments demonstrate that it is only correlated with temperature, as in next Section. 

Most of currently existing tremor estimation algorithms have been developed for cancelling physio- 
logical tremor in human-machine interfaces. Physiological tremor has another aetiology than patholog- 
ical tremors, and possesses also different characteristics: it is barely visible to the unaided eye, and has 
frequency between 8 and 12 Hz [25]. Physiological tremor is only symptomatic during high precision 
tasks, thus it is typically employed during hand held surgery. According to this, successful algorithms 
for cancelling physiological tremors such as the classical Weighted frequency Fourier Linear Combiner 
(WFLC) [16], Bandlimited Multiple Fourier Linear Combiner (BMFLC) [17], or Double adaptive BM- 
FLC [26], may not constitute the optimal solution for estimation of pathological tremor parameters. In 
fact, in [9], Rocon and colleagues employ a preliminary filtering stage that eliminated volitional motion 
from the input signal before feeding it into a WFLC, which was employed to track pathological tremor. 
A similar approach has been recently described in [27]. 

In this paper, we present a two stage algorithm that estimates instantaneous tremor parameters from 
wearable gyroscope data. As above mentioned, information about current tremor amplitude and fre- 
quency will be employed to drive a neuroprosthesis for tremor suppression based on application of se- 
lective biomechanical loads. We employ MEMS gyroscopes to record joint kinematics. Next, tremorous 
patterns are extracted from the recorded motion, based on the frequency distribution of voluntary and 
tremorous components of movement. The second stage of the algorithm estimates tremor amplitude and 
frequency from the tremorous component of joint motion. Robust, precise and low delay estimation of 
tremor parameters is achieved: average amplitude estimation error is 0.001 ± 0.002 rad/s, and frequency 
estimation agrees with spectrograms. 

This paper is organized as follows. First, Section 2 presents selection of sensors, patients, and clin- 
ical and functional tests to be recorded. Next, in Section 3, we introduce the architecture of our two 
stage algorithm, and a series of estimators that will be evaluated for voluntary movement estimation 
and tremor modelling (each of the stages). The performance of these algorithms when tracking tremor 
parameters in real data is compared in Section 4; which serves to define an optimal architecture for the 
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algorithm. Finally, we end the paper with a discussion on the results, and some conclusions summarizing 
major achievements. 

2. Methods 

In order to assess time varying tremor parameters, we record hand motion by measuring wrist flexion- 
extension. We decide to measure wrist tremor because: (1) tremors are more explicit at distal joints [1], 
(2) wrist tremor has the largest impact on disability [28], and (3) it constitutes, together with finger 
tremor, the most studied tremor in clinical literature. 

Wrist rotation is obtained with two solid state gyroscopes, placed on the distal part of the forearm, and 
on the hand. Figure 1 . Fixation on soft tissues is avoided in order to eliminate the undesired oscillations 
they create, and their intrinsic low pass filtering behavior [29]. Sensors are fastened with adhesive tape. 

Figure 1. Placement of MEMS gyroscopes (red boxes) for recording wrist flexion-extension. A differ- 
ential measurement directly provides wrist rotation. 




Proper alignment between gyroscope axes and wrist joint is ensured before recording. Wrist flexion- 
extension is simply calculated subtracting the rotation measured with the hand gyroscope from the fore- 
arm rotation. Gyroscope bias is compensated online, based on its correlation with temperature, which is 
measured by the sensor itself. Correlation between gyroscope bias and temperature in three axes of one of 
the inertial measurement units (IMUs) we employ is shown in 
Figure 2. Regression coefficients for the three gyroscopes in one of the chipsets are given in 
Equations 1 to 3, where haxis represents the bias for each axis, and T temperature in °C. 

fe^ = -0.000154 ■T + 0.00562 (1) 
by = 0.008506 ■ T - 0.35245 (2) 
6, = 0.038312 -T- 0.13358 (3) 

Sensors selected are off the shelf inertial measurement units (IMU) provided by Technaid S.L. They 
are compound of triaxial accelerometers, gyroscopes, and magnetometers. The low weight of these 
IMUs, around 40 g, makes them an optimal solution for our application, as tremor would change its 
characteristics if a larger mass were attached to the limbs [7]. Moreover, their small size (27 mm) does 
not interfere with user's movements. As above mentioned, only gyroscope recordings is considered. 
IMU sensors communicate with a hub trough Controller Area Network (CAN) bus. The host PC receives 
data from the hub through a USB interface at a sampling frequency of 1 kHz. 
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Figure 2. Relationship between gyroscope offset and temperature. Top plot shows tem- 
perature variation measured by the IMU. Bottom plot shows correlation between offset and 
temperature for X (blue), Y (green) and Z (red) axis gyroscopes. High correlation 
is observed. 
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2.1. Selection of Patients 

A group of patients affected by the most common pathologies that cause tremor was recruited for this 
study. The patients group was compound of four men and one woman, ages ranging from 48 to 74 years 
old (average 63 years old). All of them agreed to participate in the experiments, and gave their written 
informed consent. Patients had been previously diagnosed by the neurological personnel of the Hopital 
Erasme at Brussels; their clinical features are depicted in Table 1 . Patients kept their regular medications, 
to avoid tremor fluctuations related to changes in therapies. Tremor severity was assessed according to 
Faher scale [30]. The Ethical Committee of Hoptial Erasme gave ethical approval for this study. 

2.2. Clinical and Functional Tasks 

We selected four tasks that are relevant either from the clinical or usability analysis standpoint. Three 
of them are employed by neurologists to activate the different types of tremor [19], whereas the fourth 
one serves to evaluate patients' ability to perform their activities of daily living [31]. 

During the whole session, patients were sitting, with both arms comfortably resting on their lap. 
Selected tasks were: 

1. Arms outstretched: The patient is asked to elevate both arms and hold them against gravity with 
fingers abducted, hands in supination, during 30 s. This task is typically employed to activate 
postural tremor. 




0.5 1 1.5 2 2.5 3 

Time (min) 
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2. Finger to nose: The patient is asked to alternatively touch his nose and knee with the tip of his/her 
finger during 30 s. The patient must keep contact with nose and knee during a few seconds. This 
task is typically used to activate kinetic tremor. 

3. Rest: The patient is asked to keep both arms resting on the lap during 30 s. The elbow is flexed 
around 90°. This tasks is typically used to activate rest tremor. 

4. Pouring water into a glass: The patient is asked to pour 20 cl water from a standard bottle into a 
regular glass. The patient could choose how to perform the task, i.e., how to hold the bottle and 
the glass. This task is selected for functional and usability analysis. 



Table 1. Clinical features of enroled patients. 



Patient number 


Medical History 


Tremor 


Body segments affected 


Frequency 


Grade 1 


01 


Essential Tremor 


Postural 


Right upper limb 


7 Hz 


2 






Kinetic 


Right upper limb 


4 Hz 


2 


02 


Paraneoplastic Syndrome 


Kinetic 


Upper/lower limbs 


5-6 Hz 


1,2 






Postural 


Right/left inch 


2-3 Hz 


1 


03 


Idiophatic Parkinson 


Rest 


Upper limbs 


3-4 Hz 


1,3 






Postural 


Left hand 


6 Hz 


1 


04 


Extrapyramidal Syndrome 


Rest 


Upper limbs 


3-4 Hz 


2 






Postural 


Upper limbs 


3-4 Hz 


1 


05 


Essential Tremor 


Postural 


Right upper limb 


7 Hz 


2 






Kinetic 


Upper limbs 


4 Hz 


2 



^ If there are two figures the first one refers to the right limb and the second one to the left limb. If not, 
tremor severity is the same for both limbs. The grades are: 0: Absent, 1: Discrete and infrequent, not 
disturbing for the patient, 2: Mild, tremor amplitude is discrete but persistent, or tremor amplitude is 
mild but its presence is intermittent, 3: Intense, the tremor interferes in some activities, its amplitude 
is mild but it appears most of the time, 4: Severe, the tremor interferes in most of the chores, its 
amplitude is high, and it appears most of the time. 

3. Real-Time Estimation of Instantaneous Tremor Parameters 

This section presents our two stage algorithm for real-time modelling of tremor. It relies on two 
assumptions: (1) pathological tremors and voluntary motions have different frequency distributions, and 
(2) pathological tremor constitutes, from a signal processing standpoint, additive noise superimposed to 
volitional movement [32]. 
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The data we collected yields that pathological tremors occur in a frequency band higher than voluntary 
motion, Figure 3, which is in agreement with the literature. For example, one large study carried out with 
young and elderly healthy people and patients suffering from tremor, demonstrates that both groups are 
able to perform tracking tasks with a frequency up to 2 Hz, the bandwidth decreasing with age [33]. Also 
in [34], the authors perform a large spectral analysis of twenty four activities of daily living, showing that 
most of them involve wrist notion in a frequency range around 1 Hz, being the predominant frequency 
components between 0.48 and 2.47 Hz. 

Figure 3. Power spectral density of wrist rotation during a finger to nose task performed by patient 
01. It is observed that voluntary movement (below 2 Hz) has considerably more energy than tremorous 
movement (centered around 5.5 Hz). Dashed red line separates energy attributed to voluntary (left) and 
tremorous motion (right). 




Frequency (Hz) 

In addition, pathological tremors are reported to occur at higher frequencies, typically in 
the 3-12 Hz band [1,19]. Moreover, there exists a relationship between the underlying pathology and 
tremor frequency. For example. Parkinsonian tremor frequency lies within 4-7 Hz, cerebellar tremors 
manifest between 4 and 6 Hz, and essential tremors broaden to basically the whole 3-12 Hz band [1]. 

According to this, it is possible to extract voluntary and/or tremorous motion from kinematic data time 
series based solely on their different frequency contents, for example with a forth and back recursive 
digital filter that removes one of them from the original signal without causing phase distortion, i.e., 
delay. Figure 4. Although this approach is not real-time implementable, it sets the basis for our two 
stage algorithm for estimation of tremor parameters described in next section. 

However, at the moment of developing a real-time algorithm, we must take into account that power 
of tremorous component of motion is considerably smaller than that of volitional origin. This makes 
estimation algorithms based on gradient like approaches tend to converge towards the voluntary compo- 
nent, making them unsuitable for direct tremor modelling [33]. Therefore, we need to first isolate the 
tremorous component of motion. This constitutes the idea of our two stage algorithm: to first generate an 
estimation of tremorous motion (after stage 1), and next feed this signal into an adaptive algorithm that 
provides instantaneous tremor frequency and amplitude (stage 2), Figure 5. To generate an estimation 
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of tremor, we employ a tracking algorithm to estimate the voluntary component of motion. Next, based 
on the fact that tremor alters volitional motion in an additive manner [32], we remove the estimated 
voluntary movement from the input signal, obtaining an estimation of tremor. 



Figure 4. Separation of voluntary and tremorous components of movement by means of 
recursive digital filters. Top left figure shows the original signal (black) and voluntary 
movement (red) obtained with a zero phase low pass filter, fc = 2 Hz. Bottom left 
plot shows tremorous movement obtained by subtracting voluntary movement from the orig- 
inal signal. Right plots show power spectral densities of voluntary (top) and tremorous 
components (bottom). 




time (s) frequency (Hz) 



Figure 5. Block diagram of the two stage algorithm for real-time estimation of tremor parameters. First, 
we generate an estimation of the voluntary component of motion, which subtracted from the total move- 
ment yields an estimate of tremor. Afterwards, in stage two, and adaptive algorithm tracks instantaneous 
tremor parameters. 
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This section reviews a series of algorithms that are employed to develop the two stage algorithm. First, 
we present a number of techniques to track voluntary movement based on its lower frequency content. 
Next, we introduce two adaptive algorithms that track an input signal based on Fourier modelling and the 
Least Mean Square (LMS) recursion, a gradient-like approach [35], and a Kalman Filter that estimates 
pathological tremor amplitude. 
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3.1. Stage 1. Voluntary Motion Tracking 

Two types of tracking algorithms for real-time estimation of voluntary movement will be evaluated. 
As discussed above, voluntary movements are assumed to be performed between 0 and 2 Hz. Therefore, 
the tracking algorithm must be designed to neglect any component of the movement over 2 Hz. 

Voluntary movement is modelled as a first order process. If we consider a Taylor series that represents 
voluntary movement, we can neglect the second derivative x if either the sampling period Tg or the 
acceleration itself are small. Equation 4. In our case both assumptions are satisfied: the sampling period 
is 1 ms, and the maximum of x is 2.54 ■ 10~^ rad-s^^, 4 orders of magnitude smaller than average 
angular velocity. 



g-h Filters g-h filters are simple recursive filters that estimate future position and velocity of a variable 
based on first order model of the process. Measurements are used to correct these predictions, minimiz- 
ing the estimation error. Traditional applications of g-h filters are radar tracking and 
aeronautics [36]. The general form of a g-h filter is: 



Equations 5 and 6 are designated as update, tracking, or filtering equations. They estimate the current 
position and velocity of the variable, Xk,k, ik,k, based on previous predicted position and velocity, Xk,k~i, 
Xk,k-i, taking current measurement into account. Confidence on measures is weighted by gains 
and hk- Equations 7 and 8 are called prediction equations because they provide a prediction of future 
position and velocity, Xk+i,k, ik+i,k, based on first order dynamic model of the variable. As g-h trackers 
consider a constant velocity model, predicted velocity Xk+i,k is equal to the current one, Xk,k- 

g-h filters are affected by two error sources [36]: (1) the lag, dynamic, bias or systematic error, which 
is related to the constant velocity assumption, and (2) the measurement error, which is inherent to the 
sensor and measurement process. Typically, the smaller and are, the larger is the dynamic error 
and the smaller are the measurement errors [36]. Therefore, in designing a g-h tracking filter there is 
a degree of freedom in choice of the relative magnitude of the measurement and dynamic errors. To 
simplify the selection of gains, we consider two filters that are optimal in some sense. These filters are 
the Benedict-Bordner Filter and the Critically Dampened filter, described next. 



rp2 

X{t) = X{tn) + TsX{tn) + -fx{tn) + ... 



(4) 



Xk,k — Xk,k~l + QkiUk — Xk,k~l) 
ik,k = ik,k-l + T^iUk — Xk,k-l) 
■^k+l,k ■^k,k ~l~ TgXk^k 



(5) 
(6) 
(7) 
(8) 



Xk+l,k — Xk,k 



Sensors 2010, 10 



2138 



Benedict-Bordner Filter: The Benedict-Bordner Filter (BBF) minimizes the total transient error, 
defined as the weighted sum of the total transient error and the variance of prediction error due to 
measurement noise errors [37]. The BBF is the constant g-h filter that satisfies: 



h 



9 



(9) 



2-^7 

As Equation 9 relates g and h, the BBF has one degree of freedom. Because for g-h filters in- 
creasing the value of g diminishes the transient error, the larger g, the higher frequencies the BBF 
tracks. 

Critically Dampened Filter: The Critically Dampened Filter (CDF) minimizes the least squares 
fitting line of previous measurements [36], giving old data lesser significance when forming the 
total error sum. This is achieved with weight factor 6. Parameters in the g-h filter are related by: 



9 



(10) 



Selection of filter gain for the CDF is analogous to that for the BBF. 

Kalman Filter The Kalman filter (KF) is the most widespread estimation algorithm, and is employed 
in a large number of applications. We implement a KF that tracks voluntary movement modelled as a 
first order process. Equation 4. Therefore, state vector x()f:) is composed by the variable to be estimated, 
and its derivative. The problem is formulated as: 



T 

1 



1 0 



Xfc-l,fc-l 



(11) 

(12) 



Co variance matrices are defined taking into account the following considerations: 

1. Measurement noise covariance R(A;): as voluntary motion is the variable we are tracking, tremor 
is assumed to be sensor noise. The value of the measurement noise covariance is considered to be 



2^-2 



the average covariance of isolated tremor data; therefore R(A;) = cr^ = 0.0643 rad^-s 

Process noise covariance Q(A;): we hypothesize that process noise is related to voluntary motion 
changes due to tremor. A piecewise constant acceleration model is considered, [38]. This model 
assumes that voluntary movement undergoes constant and uncorrelated acceleration changes be- 
tween samples in the form of: 



Q 



S 

4 



S 

2 



(13) 



To select the variance of the random velocity component, cr^, we follow the recommendation 
in [38]: 0.5 ■ max^ < \a„\ < max^;. The second derivative of the raw recorded motion yields that 
max^ = 0.1042 rad-s"^. 
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3.2. Stage 2. Tremor Modelling 

State of the art tremor modelling algorithms rely on a time-varying Fourier series, which parameters 
are estimated recursively. Adaptation to the input signal is based on the LMS algorithm developed 
by Widrow [35]. As the LMS technique is a descend method that relies on a special estimate of the 
gradient [35], high energy voluntary motion must be removed first, to ensure proper tremor tracking. 
The first part of the two stage algorithm accomplishes this task. Figure 5. 

We evaluate the performance of two algorithms originally developed to track physiological tremor, 
and of a Kalman Filter we have developed to estimate tremor amplitude, the parameter that is more prone 
to change during the execution of a task, while tremor frequency keeps within a ±1.5 Hz interval around 
tremor frequency [39]. 

Weighted Frequency Fourier Linear Combiner The Weighted Frequency Fourier Linear Combiner 
(WFLC) is the most widespread algorithm for tremor modelling. It consists in an extension of the 
clas sical noise canceller presented in [3 5 ] , the Fourier Linear Combiner [40] , which also tracks frequency 
of the input signal based on a LMS recursion. Therefore, the WFLC adapts in real-time its amplitude, 
frequency and phase [16]: 



sin {rY!l=i^o?j , 
V 

£fc = Sfc - WjXfe - /ib (15) 



^r. = { > , \ (14) 

COS (r Y!1=i WoJ , M + 1 < r < 2M 



M 



) (16) 

r=l 

Wfc+i = Wfc + 2/ii£fcXfc (17) 

Equation 14 represents the time varying sinusoidal terms of the Fourier Series. Equation 15 defines 
the error to be minimized by the LMS recursion. Equations 16 and 17 represent frequency and amplitude 
adaptation. The WFLC has four parameters: the number of harmonics of the model, M, the amplitude 
and frequency adaptation gains, fiQ, and fii, and a bias weight /i^ that is included to compensate for 
low frequency errors [16,41]. The number of harmonics is typically fixed to 1, the other parameters are 
selected based on experimental data. 

Bandlimited Multiple Fourier Linear Combiner The Bandhmited Multiple Fourier Linear Com- 
biner (BMFLC) is a more recent algorithm derived from the FLC. It emerged to compensate for the 
limitations of the WFLC to track physiological tremor when two constituent frequencies [25] are clearly 
evident, or when frequency variations occur abruptly [17]. The BMFLC consists in a bank of FLCs 
that track the input signal based on multiple frequency components. Therefore, each FLC adapts its 
amplitude to the input signal, although its frequency remains constant. 

The performance of the BMFLC relies on the multiple fixed frequencies it can track. An interval is 
thus defined with the lower and upper frequency of the FLCs bank, t^o^ ^nd Uf. The number of FLCs in 
between is defined by parameter G. The BMFLC is formulated as follows [17]: 
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^ ^ sin{uQ+{uf-Uo)y^k), l<r<M 
^ cos (wo + (w/ - Wo) , M +l<r <2M 

ek = Sk- WjXfe - /ife (19) 
Wfc+i = Wfc + 2/i£fcXfc (20) 

Equation 18 represents the sinusoidal terms of the Fourier Series. Equation 19 defines the error to be 
minimized by the LMS recursion. Equation 20 represents amplitude adaptation. The BMFLC has six 
parameters: the number of harmonics of each FLC model, M, amplitude adaptation gain, /i, the lower 
and upper frequency of the FLC bank, wq, and ojf, and the number of filters in between, G. A bias weight 
/ib is also included to compensate for low frequency errors [16]. 

Although the BMFLC is not conceived as a frequency tracking algorithm, we developed an equation 
to estimate the current frequency of the input signal. Equation 21. Frequency estimation is obtained by 
weighting the contribution of each FLC to amplitude adaptation. For a first order Fourier series, it is 
expressed as: 

... = E (21) 



r= 



Kalman Filter The WFLC and BMFLC algorithms provide adaptation to the input signal based on 
special estimates of the gradient. On the contrary, Kalman Filter (KF) constitutes the optimal solution for 
estimation problems, in the sense it minimizes the covariance of a posteriori estimation error. Therefore 
the performance of WFLC and BMFLC, overall in terms of amplitude estimation, as it is the parameter 
that varies the most, can be enhanced using an adequate KF. In a similar manner to WFLC, we define 
our KF as: 



Ak,k-1 

Bk,k-i 

tThk-l 



0 
1 
0 



yfc 



_0 cos(a;o^ 
0 0 0 1 



0 
0 

1 

sin(cuo,fc,fc-i) 



lW0,fc-l,fc-l 

^fc-i,fc-i 
Bk-i,k-i 
_trk-i,k-i . 



Xfc-l,fc-l 



(22) 



(23) 



where A and B represent the amplitude of the sinusoidal terms of a first order Fourier series, and ojq the 
current tremor frequency. As tremor frequency suffers slow changes, and varies in a ±1.5 Hz interval 
around a characteristic frequency that depends on the patient, [39], WFLC frequency estimation is em- 
ployed. Therefore we define a cascade filter that consists on a WFLC that tracks tremor frequency, and 
feds it into the Kalman Filter that estimates tremor amplitude. 
Covariance matrices are adjusted as follows: 

1. Measurement noise covariance R(A;) = cr^, which has only slight impact on transient duration. 

2. Process noise covariance is defined as Q(/c) = diag(gi,i, ^2,2, 93,3, Qa^a), because state variables are 
considered to be mutually independent. 



Sensors 2010, 10 



2141 



4. Results 



This section presents evaluation of the algorithms described in previous section for both voluntary 
motion tracking, and estimation of tremor parameters. The idea is to find a unique filter setup that 
provides accurate tracking of instantaneous tremor parameters for every patient and task. To do so, first, 
we present the figure of merit that will be employed to tune each algorithm, and compare the performance 
of different candidates. Next, we summarize the results obtained with each of them. 

4. 1. Evaluation of Voluntary Movement Tracking Algorithms 

The Kinematic Tracking Error (KTE) evaluates the smoothness, response time, and execution time of 
a tracking algorithm, [41]. It is expressed mathematically as: 



where and cr^^.| are the mean and variance of the absolute estimation error, h* = \yk — x^+i^kV 
respectively. The former measures how fast the algorithm is capable of reacting when the velocity 
changes, whereas the latter quantifies the smoothness or filtering of the estimated variable [41]. Offline 
voluntary motion estimation obtained with a forth and back recursive filter is employed as the reference 
signal the estimators should track. This technique consists in filtering input data in both the forward and 
reverse directions; after filtering in the forward direction, the algorithm reverses the filtered sequence 
and runs it back through the filter, which yields precisely zero-phase distortion. 

First, we present results obtained with the optimal parameter(s) for each of the voluntary movement 
estimation algorithms presented in Section 3.1 The condition to select the optimal parameter(s) for each 
algorithm is to find the tuning that minimizes the KTE between filter estimation of voluntary motion and 
zero phase (off-line) recursive estimation of voluntary movement. Optimal filters are: 

• Benedict-Bordner Filter with g = 0.018. 

• Critically Dampened Filter with 9 = 0.990. 

• Kalman Filter with a^, = 0.0643 and al = 0.1042. 

The KTE is employed next to compare their relative performance. Table 2 summarizes KTE per kind 
of task and patient, for the optimal setup for the three algorithms. We observe that among the BBF, CDF, 
and KF, the CDF with 9 = 0.990 performs the best when tracking voluntary movement, as it provides 
the least KTE for all tasks; thus it is the approach we select for the first stage of the algorithm. 

Figure 6 shows an example of CDF and BBF estimation of voluntary movement from raw gyroscope 
recording, during an arms outstretched task performed by patient 01 . We observe that BBF estimation is 
less smooth than that of the CDF for a similar adaptation to transitory changes. 




(24) 
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Table 2. Kinematic Tracking Error (rad/s) for voluntary movement tracking algorithms organized by 
task. Table provides mean ± standard deviation. 



Algorithm 


Arms Outstretched 


Finger to nose 


Rest 


Water into a glass 


Benedict-Bordner Filter 
Critically Dampened Filter 
Kalman Filter 


0.194 ±0.058 
0.121 ±0.053 
0.169 ±0.100 


0.400 ±0.134 
0.372 ±0.118 
0.378 ±0.143 


0.147 ±0.091 
0.134 ±0.081 
0.174 ±0.129 


0.291 ± 0.083 
0.264 ±0.073 
0.312 ±0.124 



Figure 6. Comparison between CDF and BBF estimation of voluntary movement from raw gyroscope 
data during an arms outstretched test performed by patient 01. 




4. 2. Evaluation of Tremor Modelling Algorithms 

An optimal figure of merit to evaluate tremor estimation algorithms must consider the physical nature 
of the estimation error; for example if it originates from phase difference between real and estimated 
tremor, or estimation overshoots and undershoots. Traditional use of root mean square error (RMSE) 
suffers from these problems: (1) because errors due to undershoots and overshoots posses large power, 
which make them overshadow errors from interest, and (2) because the presence of delays affects the 
RMSE severely, although it does not necessarily indicate poor 
performance [42]. In this regard, the so-called filtered mean square error with delay correction, FMSEd, 
takes both phenomena into account [42], as it first aligns estimated tremor with the reference signal, and 
afterwards computes the delay corrected estimation error. The FMSEd is defined as follows: 

FMSEd = ^J{su-h_iy (25) 

where Sk represents the input tremor signal to be estimated, and stands for the delay compensated 
tremor estimation. Instantaneous delay dk is calculated offline by means of an adaptive algorithm that 
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minimizes the mean square error function based on a LMS-like recursion, [43]. As mentioned before, 
once instantaneous delay is obtained, RMSE between delay corrected estimation of tremor and the ref- 
erence signal obtained is computed, providing the FMSEd. 

First we select filter parameters so that they minimize FMSEd, and afterwards we compare the per- 
formance of each filter presented in Section 3.2 Optimal filters are: 

• Weighted Frequency Fourier Linear Combiner with fiQ = h ■ 10"'^, /ii = 2 • 10~^, fih = ^ ■ 10^^, 
M = 1, /o = 6 Hz. 

• Bandlimited Multiple Fourier Linear Combiner with /i = 4 ■ 10~^, /if, = 0, M = 1, /o = 3 Hz, 

/„ = 8 Hz, G = 4. 

• Kalman Filter with /iq = 5 ■ 10-^ /ii = 1 • IQ-^, /i^ = 1 ■ lO'^, M = 1, /o = 6 Hz, R = 0.01, 

Q = (1,1, 1,1). 

Table 3 summarizes FMSEj per task and patient for the optimal setup for the three algorithms. 
FMSEd is computed as the error between applying the different tremor modelling algorithms to the 
estimated tremor, obtained as the difference between raw gyroscope signal and CDF estimation, and the 
zero phase offline estimation of tremor. We observe that among the WFLC, BMFLC, and KF, the KF 
(with a previous WFLC stage) performs the best when estimating instantaneous tremor amplitude, as 
it provides the least FMSEj for all tasks, and also is capable of tracking tremor frequency in a robust 
fashion, as shown below. Therefore, we choose to employ a WFLC in cascade with a KF to estimate 
tremor parameters in real-time. 

Table 3. Filtered Mean Square Error with Delay correction (rad/s) for tremor estimation algorithms 
organized by task. Table provides mean ± standard deviation. 



Algorithm 


Arms Outstretched 


Finger to nose 


Rest 


Water into a glass 


Weighted Frequency PLC 
Bandlimited Multiple PLC 
Kalman Filter 


0.017 ±0.007 
0.007 ±0.008 
0.001 ±0.003 


0.052 ±0.023 
0.008 ± 0.019 
0.000 ±0.002 


0.014 ±0.006 
0.005 ±0.012 
0.001 ±0.001 


0.042 ±0.020 
0.006 ±0.013 
0.001 ±0.003 



Figure 7 shows tremor amplitude and frequency estimation obtained with the WFLC alone, and to- 
gether with the KF, during an arms outstretched task performed by patient 01 . We observe that the WFLC 
is capable of adapting to tremor frequency and reacting when changes occur (bottom plot), and that KF 
adapts faster than WFLC when tremor amplitude varies. This happens because of the optimal nature of 
the KF, which makes it adjust its gain continuously, on the contrary to the WFLC that has a fixed gain to 
descent the gradient. 
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Figure 7. Top plot: Comparison between tremor amplitude tracking with the WFLC (blue line), and a 
cascade algorithm compound by a KF preceded by a WFLC (red line). Bottom plot: frequency tracking 
with the WFLC (solid line), plotted against tremor spectrogram. Data corresponds to an arms out- 
stretched test performed by patient 01. 
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5. Discussion 

Previous sections presented a two stage algorithm for real-time estimation of tremor parameters. The 
first stage is in charge of separating voluntary and tremorous components of movement, based on the 
fact that tremor alters volitional motion in an additive manner [32] . Next, we generate an estimation 
of tremor as the difference between raw motion and this estimated voluntary movement. The second 
stage estimates instantaneous tremor amplitude and frequency from isolated tremorous movement. We 
proposed and evaluated three algorithms for voluntary movement estimation, and three algorithms for 
tremor modelling. 

Regarding voluntary movement estimation, we evaluated two g-h filters, the Benedict-Bordner (BBF) 
and Critically Dampened Filters (CDF), and a Kalman Filter (KF), for tracking volitional motion based 
on a figure of merit that accounts for tracking error and estimation smoothness, the Kinematic Tracking 
Error (KTE). This analysis yielded that g-h filters outperform the KF, because it results difficult to make 
it adapt quickly to changing voluntary movement patterns without tracking tremor, which translates into 
a larger KTE, Table 2. Comparing the performance of CDF and BBF, we observe that CDF outper- 
forms BBF because its intrinsic oscillatory nature, that makes it resonate almost in phase with tremor. 
This occurs because the CDF has two equally spaced zeroes; thus it behaves as a critically dampened 
system [36]. This also makes the CDF react faster when changes in volitional movement appear, decreas- 
ing estimation error during transitory periods. Therefore, the CDF constitutes the optimal filter during 
both steady and more dynamically complex tasks, as demonstrated by the fact that it provides the least 
KTE during both rest and finger to nose tests. Table 2. 
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Figure 8. Block diagram summarizing the two stage algorithm for real-time estimation of instantaneous 
tremor parameters. First, a Critically Dampened Filter estimates voluntary motion from raw kinematic 
data. Next, we generate an estimation of tremor by subtracting voluntary from raw movement. At the sec- 
ond stage, the Weighted Frequency Fourier Linear Combiner estimates instantaneous tremor frequency, 
and then feds it into a Kalman Filter that tracks instantaneous tremor amplitude. 
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Algorithms for estimation of instantaneous tremor frequency and amplitude are evaluated based 
on the Filtered Mean Square Error with Delay correction (FMSEd). This metric presents the advan- 
tages of accounting for errors due to undershoots and overshoots, and those originated from estimation 
delays [42]. First, we evaluated two algorithms originally devoted to physiological tremor tracking, the 
Weighted Frequency Fourier Linear Combiner (WFLC), and the Bandlimited Multiple Fourier Linear 
Combiner (BMFLC). Physiological tremor not only has different aetiology than pathological tremors, 
but also manifests differently in terms of amplitude and frequency [25,26]. Although the WFLC has 
also been successfully employed in the context of pathological tremor estimation [9], we evaluated 
a novel approach based on a cascade algorithm compound by a WFLC and a Kalman Filter (KF). 
This algorithm raised as the best solution in the sense it minimizes FMSEd, as it puts together accu- 
rate tremor frequency estimation based on a WFLC, and an independent KF for amplitude estimation. 
Table 3. In fact, FMSEd is at least five time smaller than that obtained with the BMFLC, which con- 
stitutes the second best candidate. KF tremor estimation proves to be robust and precise both during 
steady state regime and transitory periods, while BMFLC and overall WFLC lack of that accurate adap- 
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tation during transients, mainly because of the fixed gain they employ. This is demonstrated because 
performance in terms of FMSEd degrades the most during finger to nose and water glass tests. Table 3. 

Therefore, the optimal architecture for the two-stage algorithm is constituted by a Critically Damp- 
ened Filter that estimates voluntary movement, and a cascade algorithm compound by a Weighted Fre- 
quency Fourier Linear Combiner that estimates tremor frequency and feds it into a Kalman Filter that 
tracks tremor amplitude. Figure 8. This approach provides an average tremor estimation 
error 0.001 ± 0.002 rad/s, with robust frequency tracking if compared to spectrograms. 

6. Conclusions 

This paper presented a two stage algorithm for real-time estimation of instantaneous tremor parame- 
ters. At the first stage, the algorithm separates tremorous and concomitant voluntary movement based on 
their different distributions in the frequency domain. Next, estimated voluntary movement is removed 
from raw kinematic data, in order to generate an estimation of tremor. This estimation is then fed, at a 
second stage, into a Weighted Frequency Fourier Linear Combiner (WFLC) that tracks tremor frequency, 
and into a Kalman Filter that uses WFLC frequency information to estimate tremor amplitude. 

The resulting algorithm provides accurate estimation of tremor amplitude, with an average FMSEd of 
0.001 ± 0.002 rad/s, and robust frequency estimation when compared with spectrograms. Our two stage 
algorithm has been validated with patients suffering from major pathologies that cause tremor during the 
execution of both clinical and functional tasks. Real-time tremor parameters will be employed to drive a 
"soft robot" or neuroprosthesis for tremor suppression based on Functional Electrical Stimulation. 
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